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Above two dimensions, diffusion of a particle in a medium with quenched random traps is be¬ 
lieved to be well-described by the annealed continuous time random walk (CTRW). We propose an 
approximate expression for the first-passage-time (FPT) distribution in a given sample that enables 
detailed comparison of the two problems. For a system of finite size, the number and spatial arrange¬ 
ment of deep traps yield significant sample-to-sample variations in the FPT statistics. Numerical 
simulations of a quenched trap model with power-law sojourn times confirm the existence of two 
characteristic time scales and a non-self-averaging FPT distribution, as predicted by our theory. 
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I. INTRODUCTION 

The diffusive motion of macromolecules is an essen¬ 
tial part of cellular life. It controls the speed of a large 
number of cellular processes such as signalling, assem¬ 
bly of protein complexes and molecular machines, and 
exit of mRNAs from the cell nucleus [iHS]. Recent ad¬ 
vances in in vivo single-molecule imaging have greatly en¬ 
riched our knowledge of thermally driven transport in the 
heterogeneous intracellular medium. Experimental mea¬ 
surements have generally indicated a subdiffusive behav¬ 
ior of fluorescently tagged particles in the cytoplasm and 
the nucleoplasm, as well as on the plasma membrane^ 
|6] . The origin of the observed anomalous phenomenon is 
still much debated. 

Hitherto, the statistical characteristics of measured 
particle trajectories are mostly interpreted using one of 
the three theoretical models: the fractal Brownian mo¬ 
tion with temporal correlation of particle displacements 
[ZHII], the continuous time random walk (CTRW) with 
power-law sojourn times [El [13], or the obstructed diffu¬ 
sion caused by organelles in the diffusion path [El [T5] . 
While these models present different scenarios for the 
cause of subdiffusive scaling of particle displacement with 
time, they do not take into account the quenched na¬ 
ture of the intracellular medium. For diffusive transport 
across length scales larger or comparable to the size of 
chromosomes, endoplasmic reticulum[l()) and other or¬ 
ganelles, the environment is usually static during the pas¬ 
sage of a tracked particle. In such a situation, one expects 
significant sample-to-sample variations whose statistical 
mechanical characterization is challenging [Ej- 

One of the well-known models in this context is the 
quenched trap model (QTM) defined by a set of hopping 
rates out of sites i on a d-dimensional lattice [El- 
[23], as illustrated in Fig. 1. The time constants Ti are 
quenched random variables drawn from a distribution 
Ps{t). Above the critical dimension dc = 2 for return¬ 
ing walks, it is generally expected that the QTM has 
the same scaling properties as the CTRW where n is re¬ 
assigned according to the distribution Ps(t) upon each 



FIG. 1: (color online). First-passage trajectory (blue line) 
of a diffusing particle in a quenched environment with deep 
traps (red dots). 


visit to site i. The argument was formalized in a renor¬ 
malization group analysis by Machta[E]. 

Indications that the scaling properties of the QTM and 
that of the CTRW may not be identical in high dimen¬ 
sions can be found from previous studies of the first- 
passage-time (FPT) between start and target sites in a 
confined geometry. The FPT has been suggested to be 
an important characteristic for understanding reaction 
kinetics inside a cellOES]. Consider a simple example 
illustrated in Fig. 1 where the start site is at the center of 
a sphere of radius R and the target is any site on the sur¬ 
face of the sphere. For this geometry and a power-law dis¬ 
tribution Ps{t) ~ T~^~^ of trapping times t^, the mean 
FPT was shown to scale with R as rg ~ [T51 US] ■ 

In comparison, for /r < 1, the typical FPT of the cor¬ 
responding CTRW is given by Ttyp ~ while the 

mean FPT diverges. Since the mean FPT tends to be 
dominated by rare events when the distribution has a fat 
tail, the above results suggest qualitatively different tail 
of the FPT distribution in the two models. 

In this paper, we develop an analytic scheme to char¬ 
acterize the full FPT statistics in the QTM, and to inves- 
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tigate its sample-to-sample fluctuations. Simple physical 
arguments are presented to reveal the subtle differences 
in the rare-event statistics between the quenched and an¬ 
nealed systems. The two time scales, Ttyp ~ for the 
longest trapping time on a typical FPT trajectory, and 
tq ~ for the longest trapping time in the entire sys¬ 
tem, arise naturally from the discussion. Based on these 
understandings, we propose a decomposition scheme that 
expresses the FPT distribution as an ordered sum of ex¬ 
ponentials associated with deepest traps in a given sam¬ 
ple. Variations in the strength and location of these traps 
give rise to sample-to-sample fluctuations in the tail of 
the FPT distribution which we characterize analytically. 
The analytic results are shown to be in excellent agree¬ 
ment with numerical simulations of the two models. 

The paper is organized as follows. In Sec. II we define 
the Quenched Trap Model and mention a few of its basic 
properties. Section III presents our analytic calculations 
of the FPT distribution from the point of view of rare- 
event statistics. Section IV contains results from numer¬ 
ical simulations of the two models, highlighting sample- 
to-sample variations in the QTM. A brief summary is 
given in Sec. V. 

II. THE QUENCHED TRAP MODEL 

The Quenched Trap Model can be formulated as a 
Markov process for a diffusing particle whose states are 
sites on a d-dimensional hypercubic lattice with a lat¬ 
tice constant a. The transition rate from site i to a 
neighboring site j is given by Wi^j = {2d)~^T~^. The 
hopping rate ki = out of site i can be associated 
with a site energy Vi(< 0) through the Arrhenius law 
ki = Wo exp(Ui/T), where wq is the attempt rate and T 
the ambient temperature. Of special interest is when the 
site energies Vi is exponentially distributed, 

P{V) (F<0) (1) 

Consequently, the distribution of follows a power-law, 

Psir) = (r > W(7^) (2) 

Here the exponent /r = T/Tg. Below we will choose tq = 
W(j"^ as the unit of time unless otherwise specified. 

The QTM with a power-law distribution of sojourn 
time constants (|^ offers a plausible description of the 
diffusion of a macromolecule in the aqueous cellular en¬ 
vironment. On short time scales, the molecule is trapped 
in a certain volume of its size due to either nonspecific 
binding or cage effect as in colloidal glass. It has been 
argued that the “Gumbel distribution” approximated by 
0 is often encountered when multiple factors of compa¬ 
rable strength contribute to the trapping energy Vi m- 
On longer time scales, after breaking off from the trap, 
the particle diffuses normally until it falls into another 
trap. In this scenario, the path taken by the molecule is 
simply a normal random walk, although the journey time 


may have a very broad distribution. The lattice constant 
a can be taken to be the linear size of a correlated volume 
in the medium, beyond which the local potential Vi on 
the diffusing particle changes to a substantially different 
value. The time constant tq should be chosen accordingly 
so that ([^ provides an adequate description for hopping 
between neighboring sites. 

According to Eq. Q, the site energies are not bounded 
from below. However, in a given realization of the dis¬ 
order, there will be a deepest trap at energy Vq = 
mini {Vi}. Let V ~ be the volume of the system 
and N = V/a'^ be the number of independent sites. The 
cumulative distribution of Vq is given by 

Cv{Vq,N) = l-exp(-AAe'^«/^») 

= Cy(^+lnAA), (3) 

where the scaling function C{u) = 1 — exp(—e“). There¬ 
fore the mean of Vq decreases logarithmically with N but 
its variance remains constant. Consequently, the energy 
difference between the deepest and the second deepest 
trap, which can be approximated by the difference in Vq 
in two independent realizations of the disorder, is finite 
and independent of M. The corresponding trapping time 
constant 

tq{R) = To exp(-VQ/r) ~ To{R/aY/>^ (4) 
is greater than the second longest time by a factor. 


HI. ANALYTIC CALCULATION OF THE 
FIRST-PASSAGE-TIME DISTRIBUTION 

For diffusion-limited biochemical reactions inside a 
cell, one is interested in the FPT of a molecule from its 
birth place to the reaction site. In this work, we shall 
focus on the statistics of the FPT for individual cells. 
To simplify the discussion, we shall adopt the simplest 
geometry as illustrated by Fig. 1, where the molecule 
is launched from the origin at t = 0, and examine the 
distribution F(t, R) of the FPT t to an enclosing spher¬ 
ical surface of radius R. Nevertheless, our approach can 
be adapted to more general geometries as considered by 
Benichou et al. |24j using results from the lattice random 
walk pSl [^ . 


A. The Continuous Time Random Walk 

If the sojourn times of the diffusing particle at a given 
site are not distributed in a site-dependent manner but 
instead follow a common distribution such as Eq. 0, the 
corresponding stochastic process is known as the Contin¬ 
uous Time Random Walk. The CTRW can be considered 
as an annealed version of the QTM, with the important 
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difference that the statistics of the journey time is inde¬ 
pendent of the path taken. This “subordination” prop¬ 
erty allows one to write, 

^CTRw(t, R) = '^ (5) 

N 


where A]\[{R) is the probability for a lattice random 
walker to reach the boundary for the first time in N steps, 
and /Ar(t) is the probability distribution function of the 
journey time In = Under the continuum ap¬ 

proximation at large N (which plays the role of time) 
and R, the lattice random walk is described by the diffu¬ 
sion equation which, under the substitution x —>■ Rx and 
N —>• R^N, yields the scaling solution, 

An{R) = R-^AdiNR-^), ( 6 ) 

where the FPT probability density function (PDF) Ad{u) 
is peaked around rtmax — 1 and decays exponentially at 
large u. On the other hand, fnii) exhibits a fat tail when 
the sojourn times ti are distributed according to ([^. 

A particularly interesting case is /x < 1, where the 
mean sojourn time 0 = /g tPs{t)dt diverges. In such a 
situation, the sum is dominated by the largest term 
^max whose typical value grows faster than linear in N. 
The latter can be seen from the cumulative distribution 
of tniaxi 

N 

Prob(tinax >t) = 1 - H Prob(ti < t) = 1 - (1 - f'")^ 

i=l 

~ l-exp{-Nt-f^). (7) 

In Appendix A we show that, in this case, the tail of the 
distribution of becomes identical to that of fmaxj 

fNit) « P^(tmax =t)= Nfi{l - (8) 

Equations Q and show that the distributions of 
Ttyp and oi In are both peaked around rtyp ~ and 
have the expected power-law tail beyond Ttyp. Combining 
Eqs. ([^-([^, we obtain the subdiffusive scaling [TSl I5D] . 

FcTRw{t,R) ^ J 

= < 1) (9) 

Here $(z) = duuAd{u) exTp{—u/z^) increases mono- 
tonically with z and saturates to $oo at large z. Hence 
Actrw(C-R) is peaked at 

rtyp(A) ~To(i7/a)2/^. (10) 

The corresponding cumulative FPT distribution takes 

the form, 

/ OO 

dt'F(t',R)=C{tR-^^^). (11) 

Here C{z) = dy^iy->^-^<^{y). For z » 1, C[z) ~ 


B. The FPT distribution in the QTM 


Equation ([^ integrated over t is a special example of 
a general formula for the cumulative EPT distribution, 

C{t,R) = Y,WrCr{t), ( 12 ) 

r 

where the summation extends over all possible lattice 
walks F connecting the launch site to the boundary. Here 
Wr is the probability for a path F in the lattice walk, with 
the normalization Wr — 1, and Crit) the probability 
that the total passage time exceeds t. In terms of the 
sojourn time distribution Pi(t) on site i, we may write 

Cr{t)= f (13) 


where H{t) is the Heaviside step function. 

For the CTRW, the sojourn time distributions Pi{t) are 
identical for all sites. Hence C'r(t) depends only on the 
path length N, in which case paths with the same N can 
be grouped together, leading to Eq. On the other 

hand, C'r(t) for the QTM depends on the actual sites 
visited. An interesting question is whether the disordered 
averaged Cr{t) depends only on the path length Nr but 
not its spatial trajectory. If so, the ensemble-averaged 
C{t, R) can again be cast in the form of Eq. (§. 

The answer to the above question is in the affirma¬ 
tive if, on the lattice walk F, each site visited appears 
only once, i.e., there is no return to any given site on 
the path. For this group of lattice walks, the disorder- 
averaged C'r(t) is equivalent to its annealed counterpart 
whose sojourn time distribution is given by, 

r°° ndr , f* 

{R{t))= x^e-^dx. (14) 


Here and elsewhere we use (•) to denote average over the 
quenched disorder. It is easy to verify that Eq. (14) also 
has a fat tail that decays as at large t. 

The theory of lattice random walks [551 US] can be ap¬ 
plied to calculate the distribution of returns for the setup 
illustrated in Fig. 1. For d < 2, the typical number of 
returns grows with R as (logarithmic at d = 2). 

Hence a different scheme to compute the ensemble aver¬ 
age of C{t, R) is required. On the other hand, for d> 2, 
the return distribution decays exponentially. Further¬ 
more, the majority of returns to a given site take place 
in a short section of the walk, allowing for a renormaliza¬ 
tion group treatment of their effects on Cr {t) [T5| . In the 
following, however, we will take a different route to re¬ 
organize terms in Eq. (121 to compute C{t,R), focusing 
on contributions from the deepest traps. 

The calculation in Appendix A shows that Cr{t) has 
an exponential tail with a time constant rr = 
Furthermore, under Eq. (§ at /x < 1, Tr is dominated 
by the largest trapping constant r^axT = maxjgrlTi}. 
These observations suggest that C'r(t) may be replaced 
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by the cumulative sojourn time distribution at the deep¬ 
est trap on path F. More precisely, labelling each path F 
by T-„iax,r, we rewrite Eq. as, 

C{t,R)=J2c.{t,R), (15) 

i 

where the partial sum 

C.{t,R)= Y. WrCrit) (16) 

=Tniax,r 

is restricted to paths F that go through site i with the 
largest trapping constant on F. 

To proceed further, we introduce a single-trap model 
where r, = tq for all sites except at site i located 
a distance r from the center of the sphere, for which 
Ti = TQ ^ TQ{R/a)'^. In Appendix B, we present a cal¬ 
culation of the probability that the trap is visited by a 
first passage trajectory illustrated in Fig. 1, 

2 J 2-d _ f?2-d 

w{r,R)^{l-f,)- . (l«r<i?) (17) 

Here fd is the probability of return on an infinite lat¬ 
tice and Sd = 2Tr‘^^‘^IT{d12) is the surface area of unit 
sphere. For simple cubic lattice, /s ~ 0.34. Taking into 
account multiple visits[5S], the calculation yields an ex¬ 
ponential tail for the cumulative FPT distribution with 
a renormalized time constant. 


Ci{t,R\TQ,r) ~ ■u;(r,i?)exp(-edt/TQ), (18) 


where = 1 — fd is the probability of no return. Contri¬ 
butions from paths that do not go through the trap die 
out at times much shorter than tq. 


We now compare Eq. (16) with (181 when the site 


i corresponds to the deepest trap in the system, with 
Ti = Tq. In this case, the sum in Eq. (16) includes all 
paths that go through site z, as in the sin^-trap model. 
The total time spent on other sites on the path is of the 
order of rtyp(i?) in Eq. (16) and To(i?/a)^ in the single¬ 
trap model. Hence, in both cases, the total passage time 
is well approximated by the sojourn time on the deepest 
trap. This allows us to write, for the deepest trap. 


C,{t,R) ~ Ci{t,R\T„r,), 


(19) 


where is the distance of site i to the origin. 

For the second deepest trap in the system, the same 
argument leading to (19) applies except that the right- 


hand-side of the equation contains extra contributions 
from paths that visit both the deepest and the second 
deepest trap. To eliminate such terms, we need to re¬ 
place w{r,R) in Eq. (18) by w{i,j,R) that gives the 
probability of reaching the boundary via site i with an 
absorbing site at j. Continuing the procedure to the 
(n -|- l)th deepest trap, one needs to compute the visit 
probability to the site in question in the presence of n 
absorbing sites and the system boundary, which is quite 


challenging. However, using the annealed approximation 
that, in each step, the walker has a probability p = n/N 
to run into one of these sites, we obtain the survival 
probability (1 — p)^ — exp{—nN/N) in an Wstep walk. 
Applying this estimate to first passage trajectories with 
N R^, we see that Eq. ( |l^ holds approximately 
when n < Tic = /N ~ {R/a)°^, but for weaker traps, 
Ci(t,R) diminishes rapidly. The trapping time constant 
at Tie satisfies {t/tq)~^ = ndM = {R/ad, yielding a 


cut-off time constant Ttyp(i?) given by Eq. (10). 

The above discussion yields the following approximate 
expression for the probability that the FPT in a given 
sample is greater that t, 


C'qtm(C R) — 


E 

'tTi >'7'typ 


Ciit,R\n,n), (20) 


{R) 


where Ci{t, R\Ti,ri) at t > Ttyp(i?) is given by Eq. (18). 
Each term in the sum decays exponentially with a time 
constant given by the strength of the trap. The largest 
time constant is set by the deepest trap in the system. 

Given the simple mathematical form of Eq. (20), var¬ 
ious properties of the QTM can be derived analytically. 
For example, as we show in Appendix C, the mean and 
variance of C{t, R) over different disorder realizations are 
easily computed. For 2 < d < 4, the results are given by, 

{C{t,R)) c. (21) 


{[ACit,R)r) cs pi2edr-^Sf^ 


X7 


A-d 




( 22 ) 


Here "f{p,,z) = /p dxx^~^e~^ is the incomplete gamma 
function. As expected, the ensemble averaged FPT dis¬ 
tribution © is essentially the same as that of the 
CTRW. On the other hand, its relative fluctuation satis¬ 
fies the scaling. 


([AC(t,R)]^}V^ 

{C{t,R)) 


= i?-0-2)/2 


fitR-^'d 


(23) 


where /(z) ~ for z ^ 1. The normalized fluctua¬ 
tion is proportional to (I/tq)'''^^ in the intermediate time 
regime rtyp < t < tq, and becomes of order 1 or bigger 
when t exceeds tq. 


IV. SIMULATION RESULTS 

We performed extensive kinetic Monte Carlo sim¬ 
ulations of the CTRW and the QTM on the three- 
dimensional simple cubic lattice to verify the analytic 
results presented in Sec. HI. In each run, a particle is 
released from the origin at < = 0 and performs unbiased 
random walk through nearest neighbor hops. The walk is 
terminated when the particle, for the first time, reaches 
a site at a distance greater than R from the origin. The 
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FIG. 2: (color online). Cumulative distributions of the first 
passage time from simulations of the CTRW model and the 
QTM on the three-dimensional simple cubic lattice. Here R = 
7 and /r = 0.71. Dashed lines are computed using Eq. (201 
for the specific trap configuration in each of the two samples, 
respectively. 


total passage time is given by the sum of sojourn times 
at each stop along the path. In the case of the CTRW, 
the sojourn times are drawn independently from the dis¬ 
tribution ([^. For the QTM, a set of trapping time con¬ 
stants Ti are first assigned to the lattice sites. The actual 
sojourn time upon each visit to a given site follows an 
exponential distribution with the pre-assigned time con¬ 
stant. The system sizes investigated are from i? = 7 to 
15. 



FIG. 3: (color online). Simulation data of the QTM on the 
three-dimensional simple cubic lattice at /r = 0.71. (a) Scal¬ 
ing plot of the mean cumulative FPT distribution for three 
different system sizes R — 7 (black), 11 (red) and 15 (green), 
together with numerical evaluation of Eq. ( |21[ | (dashed line). 
(b) Scaling plot of the relative fluctuation of the cumulative 
FPT distribution from simulations of the QTM at three dif¬ 
ferent sizes, along with numerical evaluation of Eq. (231. 


A. The FPT distribution 

Figure 2 shows three examples for the cumulative dis¬ 
tribution function C{t,R) at R = 7 and — 0.71. Here 
Ayp = R^^'^ = 240 in units chosen. More than 10® trajec¬ 
tories are generated to obtain accurate statistics. As seen 
from the figure, for both CTRW and the QTM, C{t, R) 
begins to drop from its maximum value 1 around Hyp. 
The shape of C{t, R) from the two samples in the QTM 
is quite similar for t around Ttyp, confirming that the 
most probable value of the FPT in the QTM does not 
vary significantly from sample to sample and coincides 
with its annealed counterpart CTRW. In the tail part 
of the FPT distributions, however, the two samples in 
the QTM show progressively larger deviations from the 
power-law C(t,R) ~ R'^t~f^ (dash-dotted line) that de¬ 
scribes well the CTRW data. At very long times, C{t, R) 
from the QTM exhibits the exponential decay predicted 
by the theory presented in Sec. III. The tail part of the 
FPT distribution in each case is well-described by Eq. 

Figure 3 shows our simulation data for the disorder- 
averaged C{t,R) and its relative fluctuation from 10,000 


realizations of the QTM at the three system sizes R = 7, 
11 and 15. Also shown are numerical evaluations of the 
corresponding analytic expressions presented in Sec. III. 
Excellent data collapse upon scaling of the variables is 
seen over six decades in time. Both Eqs. (21) and (23) 
are in quantitative agreement with simulation data on the 
tail side of the EPT distribution. Eor t < rtyp ~ 

Eq. (21) yields a value — 0.66 less than 1, presum¬ 
ably due to contributions from trajectories not included 
in the sum ( [20| ). The latter is also responsible for the 
discrepancy between Eq. (23) and the simulation data 
in the short time regime, where the actual sample-to- 
sample variations of C{t, R) decrease rapidly due to self¬ 
averaging. Taken together, the simulation results conhrm 
quantitatively the decomposition scheme Eq. (20) for the 
tail of the FPT distribution in the QTM. 


B. Hit map 

We have also monitored the statistics of the end posi¬ 
tion of the first passage trajectories as a function of the 
passage time t. This defines a time-dependent hit map 
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on the absorbing surface. For the CTRW, since the tra¬ 
jectories are spatially independent of each other, the hit 
map is uniform at all times apart from statistical fluctu¬ 
ations when only a finite number of hits are registered. 
However, the QTM is expected to show a correlated hit 
density pattern which evolves over time. Sites with long 
sojourn time constants cause delay to trajectories going 
through them, thereby casting their shadows on the ab¬ 
sorbing surface at short times. At very long times, how¬ 
ever, all lattice random walks will have sufficient time to 
complete their journey, and uniformity is restored. 

To construct the spatially-resolved landing probability 
on the absorbing surface, we first record H{x.\t,R), the 
number of hits registered at site x up to time t. Since 
a regular lattice is used in the simulations, different sur¬ 
face sites have a different cross section to capture the 
incoming walkers. To minimize this effect, we performed 
local averaging by collecting hits into a neighborhood in¬ 
stead of a single site on the boundary. The neighborhood 
Tx of a given boundary site x includes the site itself plus 
its nearest-neighbors and next nearest-neighbors who are 
also boundary sites. Let A^x be the number of boundary 
sites in Fx. A coarse-grained hit number is defined as, 

HcG(x|t,i?) = ^ ^ H{K'\t,R), (24) 

^ x'eFx 


Using the coarse-grained data, we further normalize 
against iJcG(x|t = oo,R) to remove the lattice effect, 
i.e., by introducing the normalized cumulative hits. 


h{x\t, R) 


HcG{^\t,R) 

HcG(x|oo,i?)' 


(25) 


Figure shows the relative spatial fluctuation of 
h{x\t, R) in a system of size i? = 15 by taking statis¬ 
tics over 10^ launches. Here 


ha{t,R) = (/i(x|t,i?))x, (26) 

Ah(t,i?) = [{h\x\t,R))^-hl{t,R)]^/\ (27) 


where (•)x denotes average over all boundary sites. For 
both the QTM and the CTRW, the relative fluctuation 
decreases as the number of hits accumulate. It is also ev¬ 
ident that the spatial inhomogeneity is much stronger in 
the QTM than in the CTRW. Ideally, one expects the hit 
pattern to saturate around tq when most of the walkers 
have completed their journey to the absorbing surface. 
Equation (251 somewhat under-estimates this terminal 


fluctuation by adopting i7cG(x|oo, i?) as the normaliza¬ 
tion. The latter contains both the lattice effect as well as 
residual fluctuations from either the quenched disorder 
or a finite number of launches. 

Figure]^ a) shows the actual hit map 


u m_ h{^\t:R) 

^ (x|Ci?)= , . I, p. 


(28) 


for the QTM sample at t = 6200 ~ 3rtyp, where 50% of 
the launched particles have arrived. At this time, max¬ 
imum density variation on the surface reaches close to 



FIG. 4: Relative hit density fluctuation in the CTRW and 
QTM. Here i? = 15 and p = 0.71. 
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FIG. 5: (color online). Earth plot of the normalized hit map 
on the absorbing surface at R = 15. Here fj, = 0.71. (a) 
one sample of the QTM at t = 6200. (b) same sample at 
t = 20000. (c) CTRW at t = 6200. (d) CTRW at t = 20000. 


50%. Figure [^b) shows the hit map of the same QTM 
sample at a much later time t = 20000 ~ lOrtyp, where 
about 80% of the particles have arrived. Although the 
amplitude of the density variation has decreased from 
Fig. [^a), the pattern of high and low densities resemble 
each other. Figures [^c) and (d) show the correspond¬ 
ing maps under the CTRW dynamics, where the density 
fluctuations are much weaker. 


V. SUMMARY AND CONCLUSIONS 

In this paper, we have shown that the QTM above its 
critical dimension dc = 2 still differs from its annealed 
counterpart CTRW in the first-passage-time statistics. 
For power-law distributed sojourn times at /i < 1 where 
the particle motion becomes subdiffusive, the tail of the 
FPT distribution for a given sample is well-approximated 
by a sum of exponentials whose time constants are asso¬ 
ciated with the deepest traps in the system. Intuitively, 
this behavior can be understood from rare-event statis- 
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tics where the passage time of each trajectory is domi¬ 
nated by the strongest trap visited by the walker. By 
grouping trajectories that go through the same trap, one 
obtains individual terms in the FPT distribution, with 
weights and a renormalized time constant that can be 
calculated by applying the theory of lattice random walk. 
The largest time constant in a given sample of linear size 

d 

R scales with i? as tq ~ 7?^*. The cutoff time constant 
Ttyp ~ for traps that contribute to the summation 
is the typical FPT in such a system. 

Detailed comparison of our analytic expressions 
against simulations of the QTM shows quantitative 
agreement in the tail of the FPT distribution. For the 
peak part of the FPT distribution, however, our treat¬ 
ment is not sufficient to produce quantitatively accurate 
results, though the predicted scaling with the system size 
R is well-obeyed by simulation data. In the future one 
may consider improvements of the single-trap approxi¬ 
mation by, e.g., replacing contributions from less strong 
traps by a annealed average with a suitable cutoff, as in 
the solution of the Random Energy Model [3T]. 

The QTM differs from the CTRW also in the spatially- 
resolved arrival probability on the absorbing boundary, 
whose relative variation can be as big as 50%. It would 
be interesting to see if such hit maps can be measured in 
experiments on, say the exit statistics of mRNAs through 
nuclear pores. With sufficient statistics, the exit pattern 
may enable detection of large-scale movement of chromo¬ 
somes inside the cell nucleus. 

The work was supported in part by the NSFC under 
grant U1430237, and by the Research Grants Council of 
the Hong Kong SAR under grant 201910. 


Appendix A: Rare-event-dominated passage time 

To justify the replacement of the total elapsed time of 
a given passage by the largest sojourn time on the path 
in our treatment of the CTRW, let us compare the distri¬ 
bution of the sum Cum = Sill that of the largest 

term in the sum, tmax = max^lC}. Denoting by P{t) 
and /Ar(t) the distributions of C and tsum, respectively, 
we have. 


N 

fN{t) = (Al) 

i=l 


where angular brackets denote averaging over the distri¬ 
bution of the ti’s. Their Laplace transforms are given 
respectively by. 


P{s) 

/v(s) 


pOO 

/ P{t)e-^*dt = (e-"‘), (A2) 

Jo 

0 ; _1 


Consider now the power-law distribution Eq. (2) at 


Wo = 1- For /r < 1, its Laplace transform is given by, 

= -. (A4) 

ni n — a 

n—l 

For Ns^ <C 1, we have P(s) ~ 1—r(l —and fnis) — 
1 — F(1 — fi)Ns^. Hence we expect/v(t) ~ for 

t > in agreement with Eq. 

In the QTM, the sojourn times on the path each 
satisfies its own distribution 


Pipi) = T^ exp(-t/Ti) 


(A5) 


with a site-dependent time constant r^, the Laplace trans¬ 
form (A3) is modified to, 


/r(s)= n 


1 




1 -I- TiS 


(A6) 


For small s, we have /r(s) ~ 1 — rps where Tp = H- 
Hence /r(Q decays exponentially at large t with a time 
constant rp. In the case when r^’s are distributed accord¬ 
ing to Eq. ([^, Tp can be approximated by the largest 
time constant Tmax.r on the path. Multiple visits to the 
deepest trap can also be treated as in Appendix B. 


Appendix B: The single trap model 

In the Main Text we introduced a single trap model 
where the hopping rate ki = Tq^ except at the trap site 
located at a distance r from the center of a sphere of ra¬ 
dius i?, where h = Tq^. We shall assume tq To{R/ap 
so that the tail of the FPT distribution is dominated by 
trajectories that make at least one visit to the trap. Here 
a is the typical distance travelled by the walker in a time 
To when outside the trap. 

Let us first revisit the problem of discrete time random 
walk on a d-dimensional hypercubic lattice, paying spe¬ 
cial attention to the presence of an absorbing boundary 
S. In each step, the walker moves from its current site 
to one of its neighbors with equal probability. The prob¬ 
ability P(x,n) that the walker is at site x after n steps 
satisfies the lattice diffusion equation, 

P(x,n-tl) = ^ P{yJ,n). (Bl) 

x'Gn.n. of X 


Introduce 


OO 

n—0 

For a walker launched from the origin, summing both 
sides of Eq. ( |BI[ ) over n yields, 

^ ^ $(x0 - $(x) =(B3) 

x'Gn.n. of X 
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Here (5x,o = 1 if x is at the launch site and zero otherwise. 
Equation (B3) is a lattice version of the Poisson equation 


whose continuum limit takes the form, 


(B4) 


The solution of Eq. (B4) with the boundary $(|x| = 
i?) = 0 is given by, 


$(x) = 


2d 1 


(lxK-2 E‘i- 2 ) 


(B5) 


Here Sd = 27r‘^/Vr(d/2) is the surface area of a unit 
sphere in d dimensions. 

We now compute the probability Q(x,n) that the 
walker, launched from the origin, arrives at site x for 
the first time in n steps. Quite generally, we may write, 

n 

P{x,n) = Q{x,n')U{x,n — n'). (n > 0) (B 6 ) 

n' —1 

Here (7 (x, n) is the probability that a walker launched 
from X returns to the same site in n steps. Summing 
both sides of Eq. (B 6 I over n, we obtain. 


®(x) =5x,o + ^(x)^C/(x,n), 


n=0 


where 




(B7) 


(B 8 ) 


is the probability that the walker visits site x before being 
absorbed by the boundary S. Applying Eq. (B7) to the 
site X = 0, we obtain the probability of return, 

1 


4'(0) = 1 - 


m' 


(B9) 


For d > 2, Pq = lim^j-^oo ‘f’(O) is finite. Hence the 
Polya’s random walk constant fd = d>(0) = 1 — I/Pq is 
less than one. Equivalently, the probability of no return 
ed = 1 — /d > 0. In general, the probability of precisely 
k returns is given by. 


fd,k = (/d)^ed. 


(BIO) 


When the absorbing boundary S' is at a finite distance, 
the above results acquire a correction which essentially 
goes down as ^ 2 -d -^here ^ is the distance to the nearest 
point on the boundary. Ignoring such corrections, we 
may approximate U{x,n) by P{0,n). Consequently, Eq. 
(|B7| yields. 




Here 


wir, R) = {1- fd) 


(x^O) 


I 


2d 1 


Sdd-2\r<^-‘^ 


(Bll) 


(B12) 


is the probability to visit a site at distance r from the 
launch site, in the presence of the absorbing sphere of 
radius R. G oing back to the lattice walk, the divergence 
of Eq. (B12) at small distances should be truncated when 


r becomes comparable to the lattice constant. 

We now return to the single-trap model where the so¬ 
journ time t spent by the walker upon each visit to the 
trap satisfies the distribution: 


Pit)=TQ exp{-t/TQ). 


(B13) 


Under the assumption that tq is much greater than the 
typical FPT to the boundary when trap is not visited, we 
may write the FPT distribution of the single-trap model 


as. 


00 K 

Fi{t) ~ [1 - w{r, R)]S{t) + '^Wk{r, R)(^s(j: - ■ 






(B14) 

Here the angular bracket indicates average over the so¬ 
journ times {tj} distributed according to (B13), and 


Wk(r,R) = {fd) ^edw{r,R), {k > 0) (B15) 

is the probability that the trap is visited exactly k times 
before the walker reaches the boundary S. 


The Laplace transform of Eq. (B14) is given by, 

00 k 

Fi (s) = 1 - w{r, m + Y. ^k{r, R) J] )• (B16) 

k=l j=l 


From Eq. (B13) we obtain (e = 1/(1 -I- stq). With 
the help of Eq. (B15) we then have, 


Fi{s) = I — w{r, R) + edw{r, R) -. (B17) 

ed + STQ 


Consequently, 


Fi{t) ~ [1 - w(r,R)]5{t) F w{r, R) — exp{-edt /tq). 

'^Q 

(B18) 

See also Ref. [23] for a general discussion on the de¬ 
composition. Integrating the tail of the distribution, we 
obtain Eq. ([T^. 


Appendix C: Statistics of C(t,R) in the QTM 


The sum in Eq. (20) is restricted to sites with > 


Ttyp(i?). Alternatively, we may extend the sum to all 
sites inside the sphere, where a given site i contributes a 
term 


_ j w{n,R)e if n > Ttyp(i?); , . 


0 , 


otherwise. 
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The mean and variance of Ci are given by, 


= f dT^^w{n,R)e 

= (C2) 


{[Ac,it)f) = w^{r,,R) 




7 -M+I 


“^typiR) 

= t~^w‘^{ri,R) 

.2 


V ’ Ttvnfi?)/ 


L(2e,)M'V-’n,p(i?). 


^ V'"’rtyp(i?) 


(m, 


edt 


(C3) 


Here "^{p^z) = dxx^ “ is the incomplete gamma 
function. Since the r^’s are independently chosen from 


the distribution Eq. (2), the Cj’s are also statistically 
independent. Hence the mean and variance of C{t, R) = 
J2i Ci{t) are given by, 


(C'(t,i?))~ [ (f-xt '"ui(|x|,i?)4:7( 

J\yi\<R ^d ^ 


edt \ 

TtypiR))' 

(C4) 


J |x|<i? 


> |x|<i? 

{‘2ed)‘ 




r7 M 


'Typ(^) 


(C5) 


where we have replaced summation over i by integral over 
space inside the sphere. For 2 < d < 4, carrying out the 
integrals over x yield Eq s. pTj ) and (22|. For d > 4, 
the integral over x in (C51 is dominated by contributions 


close to the origin, and hence the lattice cut-off should 
be considered. 
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